Molecular hydrogen in the cosmic recombination epoch 
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The advent of precise measurements of the cosmic microwave background (CMB) anisotropics 
has motivated correspondingly precise calculations of the cosmic recombination history. Cosmic 
recombination proceeds far out of equilibrium because of a "bottleneck" at the n = 2 level of 
hydrogen: atoms can only reach the ground state via slow processes: two-photon decay or Lyman-Q 
resonance escape. However, even a small primordial abundance of molecules could have a large effect 
on the interline opacity in the recombination epoch and lead to an additional route for hydrogen 
recombination. Therefore, this paper computes the abundance of the H2 molecule during the cosmic 
recombination epoch. Hydrogen molecules in the ground electronic levels X^E^ can either form from 
the excited H2 electronic levels B^E^ and C^IIii or through the charged particles Hj, HeH^ and 
H~. We follow the transitions among all of these species, resolving the rotational and vibrational 
sub-levels. Since the energies of the X^S^-B^Ej (Lyman band) and X^E^-C^H„ (Werner band) 
transitions are near the Lyman-Q energy, the distortion of the CMB spectrum caused by escaped 
H Lyman-line photons accelerates both the formation and the destruction of H2 due to this channel 
relative to the thermal rates. This causes the populations of H2 molecules in X^E^ energy levels 
to deviate from their thermal equilibrium abundances. We find that the resulting H2 abundance is 
10~^^ at z = 1200 and 10"^"^ at z — 800, which is too small to have any significant infiuence on the 
recombination history. 

PACS numbers: 98.70.Vc, 95. 30. Ft, 98.62.Ra 



I. INTRODUCTION 



The era of percent-level precision cosmology started 
with the exquisite measurements of the cosmic microwave 
background (CMB) anisotropics by the Wilkinson Mi- 
crowave Anisotropy Probe (WMAP) satellite [H. The 
CMB anisotropics are a very useful tool for cosmologists 
for two reasons. First, the shapes and normalizations 
of the temperature and polarization spectra are sensitive 
to a host of cosmological parameters 0]. Second, the 
physics underlying the CMB power spectrum is thought 
to be well understood. It can be calculated by linear 
perturbation theory of the Einstein and Boltzmann equa- 
tions around a homogeneous, isotropic background 0t3|; 
the perturbation equations can be solved rapidly by mod- 
ern numerical codes that have achieved agreement at the 
0.1% level in code comparisons However, to solve 
these equations one needs to know the number density of 
free electrons as a function of redshift ne{z), the so called 
recombination history, which enters into these equations 
through the Thompson scattering of photons from free 
electrons. 

The first cosmological recombination calculations were 
carried out more than 40 years ago 0, [1], showing the 
importance of non-equilibrium hydrogen recombination 
because of the high optical depth of the Lyman series 
lines in the early universe. A hydrogen atom can only 
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reach its ground state from the H i 2s level or 2p levels 
via two-photon decay and redshifting out of the Lyman-a 
line, respectively. The early analyses assumed Boltzmann 
equilibrium of all n > 2 levels of hydrogen, and thus 
had to follow only ionized hydrogen H"*" -I- e~, excited 
hydrogen H*(n > 2), and ground-state hydrogen H(ls). 

To obtain He (z) to high accuracy, it is necessary to in- 
clude additional physics. Thus theorists have considered 
helium recombination deviations from Boltzmann 

equilibrium for the n > 2 levels of h ydroge n p^| - [23j : a 
host of two-photon processes [T^. [isl. |23 - I29j : the trans- 
port of photons near Lyman-a due to multiple resonant 
scattering [29l - [37t ; and cross-talk amon g va rious lines and 
the photoionization continuum IssI |39| . 

The workhorse recombination code Recfast, used for 
WMAP parameter constraints, was a fitting function 
to such non-equilibrium calculations including all of the 
physical processes recognized as important in the year 
'-^2000 [40I IZlj . and there have been some subsequent 
updates [i^l- Recfast was sufficiently accurate for the 
observations of its time, however to fully take advan- 
tage of the power of the Planck satellite data (launched 
2009) it is important to find ne{z) to the sub-percent 
level [4^ 113. This realization triggered a fiurry of pa- 
pers considering a host of new phenomena that could 
affect the recombination history to the percent and sub- 
percent level, culminating in two new publicly available 
codes that properly treat the radiative transfer effects in 
hydrogen and helium recombination [4^ |45| . 

The current paper is a continuation of the same ef- 
fort to reach the required level of accuracy. We consider 
how the formation and destruction of hydrogen molecules 
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(H2) in the X^E+, B^S+ and C"'^n„ electronic states can 
change the recombination history. The reason that H2 
might be able to change the recombination history is 
that the Lyman and Werner bands (X^I]+-B'^E+ and 
X^I]+-C^n„) are near the Lyman-a energy (/ifLya — 
10.2 eV). Thus the excitation, de-excitation, photodis- 
sociation, and photoassociation of the H2 molecule can 
shuffle photons between the red and the blue sides of 
the Lyman-a line. In an expanding Universe, a pho- 
ton redder than Lyman-a is likely to simply redshift and 
eventually become a part of the far-infrared background, 
whereas a photon bluer than Lyman-a will redshift into 
the Lyman-a frequency and excite a ground-state hydro- 
gen atom (which at 2 > 900 would have been likely to 
be photoionized). At an order of magnitude level, one 
would expect this H2-mediated redistribution to become 
possibly significant if the net optical depth in the Ly- 
man and Werner bands tlw were of order 10""^, which 
for the recombination-era density and Hubble rate, and 
total oscillator strengths of order unity, would require an 
abundance x[H2] = n(Il2)/nH ~ 10""'^^. Therefore we 
are interested in even tiny quantities of molecular hy- 
drogen. It is worth noting that some models of early 
Universe chemistry have found a:[H2] > 10~^^ [46l . |47| 
during the recombination epoch with simplified (i.e. not 
level-resolved) reaction networks. 

The calculation of the abundance of hydrogen 
molecules has already been considered by many authors 
but for a different cosmological goal, that is to assess the 
effect of the II2 molecule on the cooling of metal-free ga s 
and its implications for primordial star formation j48l - l56l | . 
Since direct radiative association to the X^E+ electronic 
level, i.e. 2H— 7>H2(X-'^E+) -I- 7 is forbidden, two sepa- 
rate hydrogen atoms must reach the ground state of H2 
through an intermediate route. For the case of the post- 
recombination era when fcBlcMB < 0.2 eV the accessible 
routes are through the 48] and (49,, 50^ intermedi- 
ate states. Indeed, complex reaction networks have been 
constructed to follow hydrogen chemistry [H, [H, IstI - 
[60j . These have identified in particular the significance 
of the recombination-induced CMB spectral distortion 
[17I IgU Ig^ in controlling pregalactic photochemistry, the 
importance of rate coefficients [63|-|65|, and the impor- 
tance of following transitions among the various rota- 
tional and vibrational levels of the ion at z < 500 

However, at the redshift of interest for this paper, 
z ~ 1000, there is another route for the formation of 
hydrogen molecules: the inverse Solomon process [66| . 
At this era there are enough ultraviolet photons (both 
blackbody photons and spectral distortion photons) to 
facilitate the photo-attachment of two hydrogen atoms 
into an excited H2 molecule in one of the rovibrational 
levels of either the B^E+ (Lyman band) or C^H„ (Werner 
band) electronic states with energies ^ 10 eV. The ex- 
cited H2 molecule will re-emit the photon and decay to 
either a bound H2(X"'^E+) molecule, or to the continuum 



of the X level (i.e. to two H atoms). In equation form, 
2H + 7 o H2(BiE+ C'Uu) o H2(Xil]+) + 7. (1) 

This mechanism and the charged-particle processes (H~ , 
HJ, and HeH+) control the H2 abundance at high red- 
shift. The possible effect on hydrogen atom recombi- 
nation is the main focus of this paper. We note that 
Ref. [66|] found only a small production of H2 via this 
mechanism, but they did not include the spectral distor- 
tion photons in their rate coefflcient and hence the total 
rate of H2 production could be many orders of magnitude 
larger. Of course, the same spectral distortion also drives 
H2 photodissociation ~ the left arrows in Eq. ([T]) - so the 
net effect on the H2 abundance requires a detailed calcu- 
lation. Deviations from thermal equilibrium abundances 
arise not from the amplitude of the ultraviolet photon 
spectrum, but the way in which its peculiar shape beats 
against the forest of H2 lines and dissociation continua. 
Since we work at z > 800 we will not distinguish the 
matter versus radiation temperature in this paper. 

This paper is organized as follows: in Sec. |TT] we write 
down the rate equations for the bound-bound and bound- 
free transitions. These equations are then solved in the 
steady state approximation and the results are presented 
in Sec. IIIII We discuss the size of the H2 abundances 
found in the previous chapter on the absorption of Ly-a 
photons and conclude in Sec. IIVI 

II. ABUNDANCE CALCULATION 

In the standard hydrogen recombination calculation 
[4lj one follows the evolution of several hundred energy 
levels of the Hi, Hei, and Hen atoms by including all 
bound-bound and bound-free transitions and treating the 
radiative transfer of line photons in the expanding uni- 
verse using the Sobolev approximation [ssl l67j . A similar 
treatment can be used for the H2 molecule. 

This section is organized as follows: we begin with a 
description of the reactions included (Sec. Ill A[) and then 
turn to the rate equations f Sec. Ill B"|) and the steady-state 
approximation (Sec. Ill C|) . Finally, we describe our model 
for the molecular data (Sec. Ill Dp and the radiation field 
(Sec. Ill E[) . which is required in order to evaluate the 
transition rates among H2 levels. The inclusion of the 
charged-particle reactions is described in Sec. IIIFI 

A. Basic reactions 

Here, we consider the ground X and the excited B and 
C electronic states of the hydrogen molecule. The latter 
are technically not bound since they can undergo spon- 
taneous radiative dissociation, but they are long-lived. 
We designate them with their rotational (J) and vibra- 
tional (v) quantum numbers; in the case of the C levels, 
which are A-doubled, it is necessary to describe the par- 
ity as either vector-like [0+, P = (— 1)"^] or axial [C~, 
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P = — (— l)"^]. We consider only the bound energy lev- 
els of these states up to rotational quantum number of 
J = 20. 

The bound-free radiative reactions involving these lev- 
els (we do not consider the collisional reactions in this 
paper) are the dipole-allowed transitions, 



H2(Bil]J) o H(ls) + H(ls) + 7, 



and 



H2(C± in„) o H(ls) + H(ls) + 7. 



(2) 



(3) 



The bound-bound reactions are the dipole-allowed tran- 
sitions 



H2(BiS+)oH2(XiS+) 



■7 



and 



H2(C±in„) f^H2(XiS+)-|-7, 



and the quadrupole-allowed transition, 

H2(XiS]+ f^H2(XlS+,l/'J') + 7• 



(4) 



(5) 



(6) 



There are in principle other quadrupole-allowed transi- 
tions; however those involving the B and C electronic 
states will be small compared to the dipole-allowed tran- 
sitions (B,C— >-X). We consider quadrupole transitions 
among the levels of the X electronic state with differ- 
ent rovibrational quantum numbers vJ because there are 
no allowed dipole decays from these levels, and hence 
quadrupole decay might be significant in comparison 
with excitation by (rare) ultraviolet photons. 

We will denote the "thermal abundance" of an H2 level 
by its abundance if the reaction 2H-H-H2 were in equilib- 
rium, i.e. 



(2Jj + l).ge.gnuc / 27rfi.2 m[H2] 



9l 

[Ei~2E 



3/2 



n(i!>)]/kT 2 



(7) 



where g^uc is the nuclear degeneracy, = 1 is the elec- 
tronic degeneracy, and = 4 is the degeneracy of an 
H atom. Note that since the early Universe is not in 
ionization (Saha) equilibrium, the choice of 2H0H2 as 
a reference reaction to define the thermal abundance is 
merely for convenience, and that achieving the thermal 
abundance does not imply full thermodynamic equilib- 
rium. 



B. Rate equations 

The rate equations for the above reactions can be writ- 
ten as: 



(8) 



where i denotes the level under consideration (we always 
resolve rotational and vibrational quantum numbers of 
H2), the overdot ' denotes a derivative with respect to 
proper time, and Xi = ni/nn where nu is the proper 
density of hydrogen nuclei (in any form - H 11, H i, or 
H2). 

The bound-bound term is given by: 



Xi\bb=- 



J>1 



(9) 



The first sum in the right hand side shows the rate of 
decrease of Xi by radiative decays to lower levels (spon- 
taneous -|- stimulated) and increase of Xi by radiative 
absorption. The second sum is for radiative decays from 
the higher levels to the level i and their inverse processes. 
Here fivji^) — f{vij + e) is the photon phase space den- 
sity on the blue side of the line. Also, Pij is the Sobolev 
escape probability [HI] (see Ref. [iil for a short deriva- 
tion) which is the probability that a photon emitted in 
the line to escape out of it via redshifting before being 
reabsorbed by a molecule in a lower j level. It is given 
by: 



P — 



1 



where the optical depth is 



■Ai 



93 



(10) 



(11) 



As a first step, we assume the H2 lines are optically thin, 
that is \Tij\ ^ 1, and therefore all the emitted photons 
will escape out of the resonance {Pij ~ 1). This as- 
sumption must be checked at the end of the calculation 
for self-consistency; later we will find it to be extremely 
good for all lines. 

Similarly, the bound-free term can be written as 



x^\bf^-x., I a,{Ef)[l + f{E,-Ef)]dEf 

/ ^{Ef)f{E, - Ef) dEf. (12) 

The first term on the right hand side is for the disso- 
ciation of the hydrogen molecule into two H(ls) atoms 
via radiative decay to an unbound vibrational state, 
and the second term is for the inverse process. Here 
Ehcc = — 1 Hartree is the energy of two separated H(ls) 
atoms with no relative kinetic energy. The functions 
ai{E)dE = Ai^f,:cc(E)dE are the Einstein coefficients 
for the decay from a bound state i to the continuum 
with energy E. The radiative absorption coefficients /?; 
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can then be calculated using the principle of detailed bal- 
ance: 



l + fUE,-E) ^ 
^HxjjJ^^ fth{Ei-E) 



.{E). (13) 



Using equilibrium thermodynamics to find the thermal 
abundance ratio (Eq. [7|) and plugging in the blackbody 

spectrum fth{E) = l/{e^/^^^ - 1) we finally find 



(2 J, + l)gegnuc f2Trh^ m[ll2] 
li,[E) = -2 ^-^—^ 



xe 



9l 

-(E,-E{,,,)/kBT 



a,[E). 



3/2 



(14) 



Here, = 4 is the degeneracy of the ground energy 
level of the hydrogen atom, is the electron spin degen- 
eracy of the bound state i and is equal to 1 for all the 
states considered in this paper as they are all in singlet 
electronic spin states. Finally, (/nuc is the nuclear spin 
degeneracy of the bound state i. It can be calculated by 
demanding that the total wavefunction change the sign 
under exchange of the two protons; this implies ^nuc = 1 
for the even- J energy levels of the X, B and C+ electronic 
states and for the odd- J states of C~. For the rest of the 
bound states the protons are in their triplet spin state, 

i-6 Qnuc — 3. 



C. steady state approximation 

We can rewrite the rate equations above in matrix 
form. To do that it is convenient to define the bound- 
bound transition matrix: 



A,j[l + f{v,j)] j<i 

l=j. 

In addition, we define the dissociation term 



(15) 



f ' a,{Ef)[l + fiE.,-Ef)]dEf, (16) 



and the source term 



/ " nnx'^^,^^MEf)f{E,~Ef)dEf. (17) 



Then it is straightforward to show that the Eqs. 
and (jl2p can be written in the compact form 



where the transition matrix T is defined by: 



(18) 




(19) 



If the transition times are much smaller than the age of 
the Universe, i.e. if the smallest eigenvalue Amin of T is 
^ H, we can take Xi — 0. Then the abundances can be 
found by the solution to the linear system: 



(20) 



In fact, since the steady-state H2 abundance is expo- 
nentially increasing at the end of recombination, we 
would like Amin to be larger than ^H, where C = 
(ilna;ss[Il2]/(ilna. The steady-state approximation is 
found to be valid until z « 810, i.e. during the por- 
tion of the recombination epoch most relevant to CMB 
anisotropics. We will find the smallest eigenvalue of T to 
become ~ H &t z = 750. 



D. Molecular data 

We require the bound-bound and bound-free Einstein 
coefficients Aij and a^. These must be calculated since 
no tabulations of the full radiative dissociation spectrum 
from the B and C states is available (Ref. [6§| gives in- 
tegrated radiative dissociation rates and mean photon 
energies, but not a spectrum as required here). 

For the energy levels and dipole transitions, we use the 
Born-Oppenheimer approximation. This may not be ac- 
curate for nearly degenerate vibrational levels of the B 
and C"*" states: these can mix if they have the same rota- 
tional quantum number J (C^ states cannot mix with B 
by parity conservation). In these cases, if e.g. the decay 
rate to a particular X rovibrational level is much greater 
for the pure B than the pure C+ state, then the mixing 
can enhance the decay rate from the energy eigenstate 
that is "mostly" 0+ [H, [tO]. However, in these cases we 
expect that the rates involving the B state would have a 
greater impact on the H2 abundance than the state. 
Therefore, the Born-Oppenheimer approximation is suf- 
ficient for the purpose of order-of-magnitude estimation 
of the H2 abundance. 

We have acquired the electronic ener gy l evel surfaces 
from the literature for the X fzil, B 72^ and C [zl 
states, and used them to construct Born-Oppenheimer 
wave functions. The electronic dipole matrix elements 
have also been obtained from the literature for X-B [tI] 
and X-C [t^ transitions. (The expressions for Einstein 
coefficients can be found in Appendix El) Note that for 
the radiative dissociation of the B and C electronic states, 
we consider decays to the vibrational continuum of the 
X electronic state. 

The dissociation coefficient ai{Ef) is a complicated 
function of the final energy Ef due to the existence of 
quasibound resonances of the H2 X^I]+ electronic state 
that ultimately dissociate via tunneling through the cen- 
trifugal barrier. This behavior can still be captured by 
the Born-Oppenheimer approximation as long as one uses 
sufficiently small steps in Ef] see Fig. [TJ (In principle, 
there is an additional requirement that the resonances 
not be radiatively broadened, but since they have no 
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a: 




0.001 0.01 
EjHartree] 

FIG. 1: The P-branch contribution to the spontaneous disso- 
ciation of H2(B'^Eu ) from the J — 14, v = 9 level, as a func- 
tion of the center-of-mass frame kinetic energy of the final H 
atoms. Note the quasibound resonance peak at Ek = 0.0009 
Hartree, which contains ~ 4% of the integrated rate. 



allowed decays this is not a problem for us.) In prac- 
tice, when computing the integrals in Eq. (IT^ we use 
the adaptive step size integrator odeint of Ref. [75j . 

We also consider the effect of adding the electric 
quadrupole transitions among different X^E+ states on 
the H2 abundances by using the results of Ref. ^tS] for 
the corresponding Einstein coefficients. 



proton mass and ilbfih^ = 0.022 is the baryon density 
parameter at the present time. 

The radiation field at E < Ei is obtained using the 
code of Ref. ^28*], which includes both CMB blackbody 
and spectral distortions. This code self-consistently fol- 
lows the absorption and emission of Lyman-series pho- 
tons in the Hi lines. In the case of two-photon transi- 
tions, it follows only the harder rather than the softer of 
the two photons. This leads to small errors at E < 
because the lower-energy photon in the decay 



H(2s) oH(ls) + 7 + 7 



(22) 



is ignored. However, at the redshifts of interest here 
{z > 800) this is not a significant oversight because 
the CMB blackbody is dominant at these low energies, 
and even at later times the distortion from redshifting of 
higher-energy photons (Lyman-a or the hard 2j photons) 
is more important than the soft 27 photons. The code 
also includes the spectral distortion due to two-photon 
decays from higher levels H(3s,3d), and due to the Ra- 
man process 



H(2s)+7 0H(ls)+7, 



(23) 



which results in a significant addition to the photons at 

E > Ei^ya- 

The code of Ref. [1^ does not track the extreme ultra- 
violet (EUV), defined here as photons energetic enough 
to ionize hydrogen, E > Ei. The reason is that the Uni- 
verse is optically thick to such photons: in < 10~ Hubble 
times these photons will be absorbed. There is, however, 
a constant stream of EUV photons being produced via 
direct recombinations to the ground state. 



H+ + e H(ls) +7EUV- 



(24) 



E. Radiation field 

The above equations require knowledge of the radia- 
tion field as a function of the photon energy, f{E). The 
relevant range of energies extends up to 0.54 Hartree 
(14.7 eV). This includes the range in which the spec- 
tral distortion is significant. It also extends to ener- 
gies greater than the ionization energy of Hi, Ei = 
13.6eV=0.5 Hartree. 

The CMB blackbody component is specified by the ra- 
diation temperature. This is T{z) = 2.728(1 -I- z) Kelvin, 
and we have not distinguished between the matter tem- 
perature and the radiation temperature (in our redshift 
range z > 800 these differ by < 0.1%). The number 
density of hydrogen nuclei is 



"h(z) = (1 + zy, 



(21) 



where Y 



0.24 is the mass weighted primordial helium 
abundance, pcr — 1-8788 x 10~^^/i^ g cm~'^ is the critical 
density of the universe at the present time, rrip is the 



The abundance of the se p hotons will thus rapidly reach 
its equilibrium value ^38]; in the limit of f{E) <C 1 so 
that we can neglect stimulated recombinations, we have 



fiE) = e 



_ -E/kBT_ 



i[R+]n,/n[R{ls)] 
{n[H+]ne/n[H(ls)]}saha 



(25) 



[The "temperature" in this equation (both in the expo- 
nential and the Saha abundance ratio) is technically the 
matter temperature T^^ since the photons are being pro- 
duced and destroyed by interaction with matter, but at 
high redshift we do not make this distinction.] 

The overall radiation spectrum is shown in Fig. [5] at 
z = 1142. 

In our analysis, we take the output of the standard 
calculation, namely the abundance of hydrogen atoms 
and the phase space density of photons as a function of 
redshift and energy, as the input in our calculation of 
the abundance of hydrogen molecule levels. Of course, 
if H2 abundances turn out to be high enough to make a 
considerable change in the recombination history or radi- 
ation spectrum, one must perform a more self-consistent 
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0.2 0.3 0.4 0.5 0.6 

E^[Hartree] 



FIG. 2: The radiation spectrum a.t z — 1142, sliowing both 
the CMB blackbody spectrum (straight line) as well as the 
full calculation including H spectral distortions. 




FIG. 3: The network of routes for the formation and destruc- 
tion of an H2 molecule that we consider in this paper. 



calculation in which the feedback of H2 molecules on hy- 
drogen recombination is taken into account properly by 
solving the rate equation for all of the species simultane- 
ously. But for the standard cosmology, we will see that 
this situation does not arise. 



F. The charged-particle processes 

We have investigated the processes that produce H2 
from neutral hydrogen atoms and radiation. However, 
there are other mechanisms that contribute to their for- 
mation and destruction, namely the H~ , H J , and HeH+ 
pathways (these have been found to be dominant at low 
z in previous works, e.g. Refs. [46l.[60|). 

All of the reactions that we consider in this paper, 
including these pathways, are shown in Fig. |31 They 
include the radiative attachment/detachment for 



H + e-oH-+7, (26) 

the radiative association/dissociation of 

H + H+ H+ + 7, (27) 

and that for HcH+ 

He + H+ HcH+ + 7. (28) 

There are subsequent nonradiative reactions that gen- 
erate H2 (and at sufficient temperature can destroy it by 
operating in reverse): the H~ channel 

H-+HoH2-^e- (29) 

and the HJ channel 

HJ+H0H2+H+. (30) 

The latter can be aided by proton exchange from HeH+ : 

HeH+ -^HoHe-^H+. (31) 

1. Positive channel: H2 and HeFt^ 

The positive ion channel (Hj") must be treated via 
level- resolved chemistry [4^ [53, |63|- Our solution to 
this is to extend the T-matrix to include the additional 
species. That is, we write 

(Th2,H2 Tjj^,H+ TH2,HeH+ \ 
Th^h. Th+3+ Th+3,h+ , (32) 
THeH+,H2 Tjj^fj+,HcH+ j 

and similarly the source vector s and abundance vec- 
tor X are extended. The model H2 molecule includes 
1435 levels (X, B, C+, C~) with J < 20. We follow all 
423 rovibrational levels of H^(X^S+). However, we treat 
the HeH+(X^E+) ion assuming Boltzmann distribution 
of the rovibrational levels at the radiation temperature, 
which is a good approximation since electric dipole tran- 
sitions are allowed and rapidly thermalize the level pop- 
ulations. The full T-matrix is thus 1859x1859. 

The sub-block of the T-matrix and s- vector involving 
HJ and HeH+ was computed in Ref. 'g^: it contains 
Eqs. (|27l) . (|28)) . and ([3T|). as well as contributions associ- 
ated with the electric quadrupole transitions among the 
rovibrational states of H^X^I]+). We used the code and 
rate coefficients of Ref. [6^ to generate the corresponding 
matrices and source vectors. Since none of our reactions 
directly connect HeH+ to H2 (they can only intercon- 
vert via H+), we set Th2,hcH+ = and Thoh+32 = 0- 
It therefore remains only to determine the sub-blocks 
Tjj^ H+ ^^'^ '^H+ H2' '^hich arise from Eq. (pOl) . These 
are related by the usual detailed balance relation, 

Rji_ ^ ■^■t,th ^ a;[H(ls)] gi ^(Ei+Ea,i,)-E,-E„+)/kBT 

R13 Xj^th a;[H+] 2gj 

(33) 
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where i represents any state of H2, j any state of Hj, 
£'h(is) —Eyi+ = —13.6 eV, and the factor of 2 comes from 
the degeneracy ratio of H versus H^" . The corresponding 
contributions to T can then be determined from Eq. ([TS]) . 

The rovibrational level-resolved forward reaction rates 
Rij for Eq. ((30)) are unfortunately not available in the 
literature. Vibrationally resolved quantum-mechanical 
rates have been calculated [z3l, however their calcula- 
tion did not resolve the rotational levels and did not 
well-sample the lowest energies required here. There is 
also the experimentally measured low-temperature rate 
of 6.4 X 10~^° cm^ s~^ [zl], which once again did not 
resolve rovibrational levels (see also the recent measure- 
ments of the reaction of D^-|-H [t^). The experiment of 
Ref. [78*1 did however show by isotopic substitution that 
the reaction mechanism is charge transfer, i.e. the two 
nuclei in the initial H J ion remain in the H2 molecule that 
is produced. Therefore the ortho- or para- nuclear spin 
character should be preserved in the reaction of Eq. (pO)) . 

Since the calculations of Ref. [t^] suggest that at en- 
ergies of several tenths of an eV the reaction H^(i^ = 
0)-|-H— )-H2+H+ is most likely to leave the final molecule 
in the v = 4 vibrational state, a simple prescription is 
to assume that (i) the rate coefficient for this reaction 
is 6.4 X 10^ ^'^ cm'^ s^^; (ii) the final vibrational state is 
J/ = 4; and (iii) the rotational quantum number (N for 
and J for H2) is unchanged in the collision. 



2. Negative channel: H 

Previous work '46', ^5^, 16011 has established that at 
z > 200, reaction Eq. (^5)) forces the abundance to 
its thermal equilibrium value because (i) the matter and 
radiation temperatures are equal, and (ii) photodetach- 
ment by blackbody photons is the main sink for H~ on 
account of the strong CMB field and the low binding en- 
ergy of H~, -B[H~] = 0.754 eV. The thermal abundance 
is then given by 



nirknT 



3/2 



,B[H-]/fcBT 



(34) 



The associative detachment reaction, Eq. (p9)) . can be 
incorporated by adding appropriate sources and sinks. 
For the sources, we add 



Si+ = fci7iHa;[H ]a;[H(ls)], 



(35) 



where "+ =" means that the quantity is added to s^. 
Here ki is the rate coefficient to level i of H2, which 
we obtain from the calculations of Ref. SOi] ; we use the 
T = 3000 K column of their table as it is most appro- 
priate for the recombination epoch and the temperature 
dependences are weak. Detailed balance implies a corre- 
sponding sink for H2 molecules: 



kiniix[Il-]x[R{ls)] 



2^1, th 



(36) 



(We may use Xi^th here since the negative species are 
catalysts and hence have no effect if H2 is at the ther- 
mal abundance; one may check explicitly that Tu in fact 
is proportional to nfj(is) with a coefficient that depends 
only on temperature.) 



III. RESULTS 

A. H2 abundance 

We now determine the H2 abundance by two methods: 
first, with only the Lyman and Werner band reactions 
as sources and sinks for H2; and second, including the 
charged particle reactions as well. The latter, of course, 
is our final result. 



1. H2 dipole and quadrupole transitions only 

The abundances of the H2 molecule in different rovi- 
brational states Xi can be found by using Eqs. p4H20|) . 

The total abundance of all H2 molecules irrespective 
of rovibrational state is shown in Fig. U) This abundance 
shows the expected rapid increase in the early stages of 
recombination as the temperature drops. At z < 1400, 
the photon phase space density f{E) begins to deviate 
substantially from the Planck spectrum. One can see 
that the sense of the resulting spectral distortion is a 
temporary decline in the abundance of H2 , followed by a 
rapid recovery as the spectral distortion redshifts below 
the Lyman- Werner bands. 

We show in Fig. [5]the ratio Xi/xi^th of the abundances 
of the states in all vibrational levels within the X^I]+ 
J — sequence of the hydrogen molecule to their abun- 
dances in thermal equilibrium at z = 1142. Here again 
"equilibrium" refers to the abundance Xi,th that would be 
obtained if the reaction H2 (X^ , i^, J) -o- 2H(ls) were 
in equilibrium at the actual H(ls) abundance (the defi- 
nition is important since the ionization fraction deviates 
from Saha) . The first two vibrational levels have approx- 
imately the same Xi/xt^th] then there is a sudden drop 
in this ratio going from the second to the third vibra- 
tional level, after which Xi/xi th increases steadily and 
approaches 1 for the weakly bound states. The physi- 
cal reason for this situation is that a photon in or red- 
ward of the H I Lyman-a line is capable of exciting an H2 
molecule only from the j/ > 2 vibrational levels (starting 
from v — 0,1 there is insufficient energy to reach the 
B or C electronic states). Thus the > 2 vibrational 
levels are rapidly photodissociated. The inclusion of the 
quadrupole transitions enables H2(X^S+) molecules to 
switch among the various vibrational levels and leads to 
a washing-out of the step in abundance versus v. 
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FIG. 4: The abundance of H2 molecules, a::[H2], as a function 
of redshift. Short dashed blue and dotted magenta lines show 
the cases where we do and do not include quadrupole transi- 
tions between X levels, respectively, and only the Lyman and 
Werner band reactions can create and destroy H2. The long 
dashed red curve shows the case when all the relevant reac- 
tions are included. The green solid curve shows the thermal 
abundances of Eq. ((T)). 
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FIG. 5: The ratios of the abundances of the II2 X^E^ levels 
with J = to their thermal values at z = 1142, plotted as a 
function of the vibrational quantum number v. Blue stars and 
magenta squares show the cases with and without quadrupole 
transitions between the X electronic levels, respectively. This 
plot included only the Lyman and Werner band reactions as 
sources and sinks for H2. 



2. Inclusion of charged particle reactions: H2 , HeH^ , and 

H- 

We now turn on our full reaction network. We display 
the abundances of the various levels of H2 by plotting 
the logarithmic abundance \ogiQ{xi/ gi) versus the level 
energy. This is a straight line in the case of a thermal 
distribution of levels. The actual result is shown in Fig. [6] 
for the X electronic states on the left and the B, C+ and 
on the right. We can see that the highly excited rovi- 
brational levels of the X electronic state are near thermal 
equilibrium, but the lower levels are underpopulated by 
~ 2 orders of magnitude. On the other hand the B, C+ 
and C~ states are overpopulated compared to the equi- 
librium abundance and they become more overpopulated 
relative to thermal for higher energy levels. This is a con- 
sequence of the shape of the spectral distortion. 

Thus far we have assumed the steady state approxima- 
tion, i.e. that all of the eigenvalues of the T-matrix are 
large compared to the Hubble expansion rate H. The 
minimum eigenvalue of T is shown in Fig. [7] together 
with the Hubble rate. We can see that for z > 800 
the minimum eigenvalue Amin is much larger than the 
Hubble rate. In fact, due to the rapid increase in the 
H2 abundance, we should really be comparing Amin to 
d\nx[}i2]/dt] the crossing of these occurs at z ^ 810, 
which is roughly where we expect the steady state ap- 
proximation to fail. As a consequence, we do not show 



H2 abundances at lower redshifts. 

In Fig. 0] we show the sum of the abundances of all 
the X, B and levels of the H2 molecule as a function 
of redshift (in practice, this sum is dominated by the X 
electronic state) . Short dashed blue and dotted magenta 
lines show the cases where we do and do not include 
quadrupole transitions between X levels, respectively, 
and only the Lyman and Werner band reactions can cre- 
ate and destroy H2. The long dashed red curve shows 
the case when all the relevant reactions are included. 
The green solid curve shows the thermal abundances of 
Eq. ([7]). We can see that the addition of new transi- 
tions consistently causes the non-thermal abundance of 
H2 molecules to increase, while always staying below the 
thermal abundance. 

We also show the abundances of the intermediate 
species in Fig. |51 The left panel shows the abundances of 
the charged particles and the right panel the excited elec- 
tronic states B and of the hydrogen molecule. Even 
though the abundances of the excited states are much 
smaller than those of the charged particles, they must be 
included since the transition rates between these excited 
states and the H2 ground state can be much higher than 
the transition rates connecting the charged particles to 
the ground state. 

A related possible formation channel for H2 that is 
mentioned in the literature is from the reaction of an 
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FIG. 6: The abundance of each H2 level plotted versus the level energy at z = 1142. In the left panel we show the abundances 
of X levels and in the right panel the B, C"*" and C~ levels. The solid line shows the thermal equilibrium abundances. 
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FIG. 7: The minimum eigenvalue of the rate matrix T (solid 
blue) as compared to the Hubble rate (dashed magenta). 



excited and a ground-state H atom [66[: 

H(Ls)+H(n = 2) f^H2(Xil]+)+7. (37) 

This reaction proceeds if the reactants approach each 
other in the vibrational continuum of a "'^Ej or ^n„ elec- 
tronic state (usually B or C), and in order to produce 
a bound H2 molecule the photon must have an energy 
E > E{Lya). Since the phase space density of photons is 



a steeply decreasing function of energy (including a step 
at Lyman-a), we would expect that the photodissociation 
of the levels of H2 is dominated by transitions through 
the discrete vibrational levels of B^S+ and C^IIu rather 
than the higher-energy continuum. One can also check 
the formation of H2 by this mechanism directly: using 



the rate coefficient 2.09 x 10 



-14 



(r/300 K) 



0.24 -T/37800K 



57| , we find a maximum production rate of H2 
molecules per H nucleus per Hubble time of 7 x 10~^^ 
at z = 1250. This channel is thus subdominant to the 
HeH+/Hj channel. All of these are small compared to 
the transitions through H2(B^S+,C"'^n„). 



B. Effect on recombination 

Since some of the X-B and X-C transitions have en- 
ergy above the Lyman-a energy, the absorption or emis- 
sion of a photon in this line by an H2 molecule adds or 
removes a photon that would have otherwise excited a 
hydrogen atom at a later time, and therefore it will al- 
ter the process of recombination of hydrogen atoms. To 
estimate this effect, we must estimate both the Sobolev 
optical depth of the various lines using Eq. PT|) . and the 
rate of production of photons in the X-B and X-C bands. 

In Fig. [5] we show the net Sobolev optical depth for 
all X~B and X-C lines. We can see that the optical 
depth of the sum of all the lines is much smaller than 
one, and therefore each individual line is also optically 
thin. This then justifies using Pij = 1 for the Sobolev 
escape probability throughout this paper. 

The net rate of emission of photons in the X-B and X- 
C bands at energies above the Lyman-a energy is shown 
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FIG. 8: Here we show the abundances of the intermediate species, H 
right panel, as a function of redshift during the recombination epoch. 



HJ, and HeH"*" in the left panel and H2(B,C ) in the 
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FIG. 9: Total Sobolev optical depth of all of the Lyman and 
Werner lines as a function of redshift. 



in Fig. [TUl We see that there was no absorption or emis- 
sion in this line at very early times when the universe was 
in thermal equilibrium. However, at redshifts of ^ 1500 
the deviation of the hydrogen molecule abundances from 
their thermal values, coupled with the CMB spectral dis- 
tortion, lead to a net emission of super Lya photons for 
both the Lyman and Werner bands. The Lyman band, 
however, started absorbing super Lya radiation later on. 
The net number of emitted photons turns out to be neg- 



FIG. 10: The net rate of emission (i.e. emission minus ab- 
sorption) of photons at -B > -E'(Lya) in the H2 Lyman and 
Werner bands, measured in photons per H nucleus per Hubble 
time. 



ligiblc (^10 and so it can only cause a very small 
change in the abundance of H atoms. 

There is one final possibility for reactions involving H~ 
and to affect recombination, namely through the neu- 
tralization reactions 



+ H+ 2H 



(38) 



11 



and the high cross section for Lyman and Werner band ab- 

sorption - H2 is not relevant for determination of the 
+ O 2H, (39) primordial recombination history and CMB anisotropics. 



which - in combination with Eqs. ([26| and (|27|) - lead 
to a net recombination. They are, however, negligible: 
even for the "upper limit" rate coefficient for Eq. (|38|) 
of 5 X 10~^ cm^ [s^l, we obtain a maximum recom- 
bination rate per Hubble time {xis/H) of 1.1 x 10^^^. 
Even this is too large since Eq. ([38]) usually leaves one 
of the hydrogen atoms in an excited level (principally 
n = 3) |8l| from which it has a large probability of being 
photoionized. Similarly, for Eq. (P^)) . the rate coefficient 
at recombination-epoch temperatures is of order ^ 10^^ 
cm'^ s""'^ (82| . which implies recombination rates Xis/H 
of order 10"!°. 
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IV. DISCUSSION AND CONCLUSION 

Since the Lyman and Werner band transition energies 
of H2 are near the H Lya energy, it is expected that the 
abundances of H2 energy levels deviate appreciably from 
their thermal abundances. This is because the photon 
phase space density has been distorted by the redshifted 
Lyck photons as in the standard hydrogen atom recombi- 
nation picture. However, it is not clear from the outset 
whether this distortion to the photon phase space den- 
sity increases or decreases the abundances of H2 levels 
compared to their thermal abundances, since the spec- 
tral distortion photons accelerate both the production 
and destruction of H2. To answer this question and ulti- 
mately to see to what extent the H2 molecules can affect 
the recombination history we have in this paper carried 
out a detailed calculation including all of the rovibra- 
tional levels of the H2 X^I]+, B^I]+, and C"'^n„ electronic 
states up to rotational number J = 20, together with the 
charged species relevant to the formation of hydrogen 
molecules, that is HJ, HeH+ and H~. We have calcu- 
lated the bound-bound and bound-free dipole transition 
rates for the Lyman and Werner bands of the hydro- 
gen molecule using the Born-Oppenheimer approxima- 
tion. Special care has been taken to find the resonances 
of the bound-free transitions. The rate equations con- 
necting the energy levels are then solved in the steady 
state approximation and the level abundances are found 
by a matrix inversion for each given redshift. 

The main result of our paper is that the shape of 
the CMB spectral distortion reduces the abundance of 
H2 compared to the thermal abundance, resulting in 
low H2 abundances throughout the recombination epoch; 
see Fig. |4l The inclusion of the quadrupole transitions 
among rovibrational levels of the X electronic state in- 
creases the H2 abundance, and adding the charged par- 
ticle processes increases the H2 abundances yet more, 
while remaining below the thermal abundance. We find 
a^[H2] ~ 10~^^ during most of the recombination epoch, 
rising to 10""'^^ at z = 800. We conclude that - despite 



Appendix A: Formulae for H2 dipole transitions 

This appendix summarizes the notation and formulae 
we use for the H2 molecular wave function and the dipole 
transition rates. 

We use unprimed coordinates [x, y, z) to denote a 
laboratory-fixed frame, and primed {x',y',z') to denote 
a coordinate system that rotates with the molecule: the 
z'-axis is taken to be parallel to the internuclear separa- 
tion vector R, and the x' axis is then chosen to lie in the 
zz'-plane. 

For the X^I]+ electronic state, the wavefunctions can 
be written in the Born-Oppenheimer approximation as 



X,uJM 



(r,R) = 



R 



y7M(R)x(r|R), (Al) 



where r = (ri, r2) is the 6-dimensional vector of positions 
of the two electrons, and x(r|R) is the electronic energy 
eigenstate associated with the X^E+ state with inter- 
nuclear separation R. A similar equation holds for the 
B^i;+ electronic state. For C^H„, one must account for 
the existence of two degenerate wavefunctions x±(r|R) 
with z'-component of electronic orbital angular momen- 
tum A = ±1. In this case the wave function is 



*c±,i.JM (r,R) = 



y+i(R)x+(r|R) 



0.j(i?) 
±y7^i(R)X-(r|R) 



(A2) 



where Y^j^f ^ spin-s spherical harmonic [8J, ISJ]- (Some 
references, e.g. Ref. \70\, use rotation matrices instead of 
spin-weighted spherical harmonics, but these are equiv- 
alent.) Note that there are two such wave functions (+ 
and — ) with opposite parity, which are degenerate (A- 
doubled) in our level of approximation. 

Normalization of the total wavefunction then re- 
quires \(j)^j{R)\^dR = 1. Since the C electronic 
states have parity u, i.e. x(— r|R) = — x(r|R), 
the total wavefunction satisfies: ^'^i j^jm (r, — R) = 
±(-l)-'+ivI/c±,.jM(r,R). 
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The Einstein coefficients for dipole transitions are dis- 
cussed in general textbooks, e.g. Refs. [H, [sl]. In our 
case, the C— >-X transition rates are as follows: for the 
electric dipole transition of a H2 molecule from an ex- 
cited state a — (C, Va, Ja) to a lower level b = {X, v^,, Jh), 
we have 



\Mab 



(A3) 



3ft4c3(2J, + i; 
where the matrix element is 

Ja Jb 

(A4) 

Here d — — eri — er2 is the electric dipole operator. The 
matrix element expands as 



XvbJbMb 



R V2R ^ 

y+i,^(R)D^^(R)±r,;i,JR)D^^(R)] , (A5) 

where the fixed-R electric dipole moment for the C —?' X 
transition is: 

Dg^(R) = -e J drld4x*x{r\R){ri+r2)xc,±i^m 



-DcxiR)- 



V2 



(A6) 



The values of Dcx{R) are calculated in |73|. For the 



B^-X transition one finds D^^(R) = Dbx(-R)R- We 
use the results of ^4] for it. 

By plugging Eq. (IA6|) into Eq. (jA5l) . separating the R 
integration from the R integration, writing x' and y' in 
terms of the fixed coordinated unit vectors (x,y,z) with 
coefficient written in the form of spin wighted spherical 
harmonics with degree s = 0, calculating the integrals 
of the three spin-weighted spherical harmonics in terms 
of the Wigner 2>j symbols and using the orthogonality 
relations of the ij symbols to do the sums over Ma and 
Mfc one finally finds: 



{Ja + l)K{C+UaJa,XubJb), and 

{2Ja + l)K{C-VaJa,XvbJb), (A7) 



where the change in angular momentum is denoted by 
the branch indices P {Jb — Ja = 1), Q {Jb — Ja = 0), and 
R {Jb - Ja = -1), and 



K{YaVaJa,Ybl'bJb) 



{R)DAB{R)4>Ya^ajSR)dR 



(A8) 



Similarly, for the B— transitions. 



{Ja + l)K{BvaJa,XvbJb) and 

^ JaK{BiyaJa,XvM. (A9) 
These equations agree with Ref. 70], appropriately re- 
stricted to the case of no B-C mixing. 
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